home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / gauss_cvf.pro < prev    next >
Text File  |  1997-07-08  |  2KB  |  71 lines
  1. ;$Id: gauss_cvf.pro,v 1.3 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       GAUSS_CVF
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the cutoff value (v) such that:
  11. ;                   Probability(X > v) = p
  12. ;       where X is a random variable from the standard Gaussian (Normal)
  13. ;       distribution with a mean of 0.0 and a variance of 1.0
  14. ;
  15. ; CATEGORY:
  16. ;       Statistics.
  17. ;
  18. ; CALLING SEQUENCE:
  19. ;       Result = Gauss_cvf(P)
  20. ;
  21. ; INPUTS:
  22. ;       P:    A non-negative scalar, in the interval [0.0, 1.0], of type
  23. ;             float or double that specifies the probability of occurance
  24. ;             or success.
  25. ;
  26. ; EXAMPLE:
  27. ;       Compute the cutoff value (v) such that Probability(X > v) = 0.025
  28. ;       from the standard Gaussian (Normal) distribution. The result should 
  29. ;       be 1.95997
  30. ;         result = gauss_cvf(0.025)
  31. ;
  32. ; REFERENCE:
  33. ;       ADVANCED ENGINEERING MATHEMATICS (seventh edition)
  34. ;       Erwin Kreyszig
  35. ;       ISBN 0-471-55380-8
  36. ;
  37. ; MODIFICATION HISTORY:
  38. ;       Modified by:  GGS, RSI, July 1994
  39. ;                     Minor changes to code. New documentation header.
  40. ;-
  41.  
  42. function gauss_cvf, p
  43.  
  44.   on_error, 2  ;Return to caller if error occurs.
  45.  
  46.   if p lt 0. or p gt 1. then message, $
  47.     'p must be in the interval [0.0, 1.0]'
  48.   if p eq 0 then return,  1.0e12
  49.   if p eq 1 then return, -1.0e12
  50.  
  51.   if (p gt 0.5) then begin
  52.     p = 1 - p
  53.     adjust = 1
  54.   endif else adjust = 0 
  55.  
  56.   below = 0
  57.   up = 1.0
  58.  
  59.   while gauss_pdf(up) lt 1.0 - p do begin
  60.     below = up
  61.     up = 2 * up
  62.   endwhile
  63.  
  64.   x = bisect_pdf([1.0 - p], 'gauss_pdf', up, below)
  65.   if adjust then begin
  66.     p = 1 - p
  67.     return, -x
  68.   endif else return, x
  69. end
  70.  
  71.